Analysis of environmental effects on Chinook salmon from the Willamette River basin
Mark D. Scheuerell
Fish Ecology Division, Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA USA, mark.scheuerell@noaa.gov
J. Tyrell Deweber
Oregon Cooperative Fish & Wildlife Research Unit, Department of Fisheries & Wildlife, Oregon State University, Corvallis, OR USA, Jefferson.Deweber@oregonstate.edu
Tom Friesen
Oregon Department of Fish and Wildlife, Corvallis, OR USA, tom.friesen@oregonstate.edu
This is version 0.17.09.27.
1 Requirements
All analyses require the R software (v3.3+) for data retrieval, data processing, and summarizing model results, and the JAGS software (v4.2.0) for Markov chain Monte Carlo (MCMC) simulation. Please note that some of the R code below may not work with older versions of JAGS due to some changes in the ways that arrays are handled.
We also need a few packages that are not included with the base installation of R, so we begin by installing them (if necessary) and then loading them.
if(!require("EGRET")) {
install.packages("EGRET")
library("EGRET")
}
if(!require("R2jags")) {
install.packages("R2jags")
library("R2jags")
}
if(!require("RCurl")) {
install.packages("RCurl")
library("RCurl")
}
if(!require("gsl")) {
install.packages("gsl")
library("gsl")
}
if(!require("reshape2")) {
install.packages("reshape2")
library("reshape2")
}
if(!require("viridis")) {
install.packages("viridis")
library("viridis")
}2 User inputs
We begin by specifying the names of three necessary data files that contain the following information:
- observed total number of adult spawners (escapement) by year;
- observed age composition of adult spawners by year;
- observed total harvest by year;
- metadata for the covariates.
The metadata for the covariates is a .csv file with the following columns:
spp: the species to which the covariate applieslife_stage: life stage at which the effect is thought to occurcovariate: type of covariate (e.g., flow, temperature)code: USGS code for data type (i.e., usually 5-digit integer)long_name: long name for the covariateshort_name: function name to derive the covariatelag_1: years to lag begin datelag_2: years to lag end datebegin: beginning date as 2-digit text (ie, mo-yr)end: ending date as 2-digit text (ie, mo-yr)location: description of locationgage: USGS gage numberflag: flag to include (1) or exclude (0) the covariateflow_scen: flag to include (1) or exclude (0) flow scenariogroup: integer indicator for life stage grouping
Let’s also define the following parameters, which will be referenced throughout the analysis.
n_yrs: number of calendar years of dataA: number of age classesM: number of covariates
## 1. file with escapement data
## [n_yrs x 2] matrix of obs counts; 1st col is calendar yr
fn_esc <- "chin_esc.csv"
## 2. file with age comp data
## [n_yrs x (1+A)]; 1st col is calendar yr
fn_age <- "chin_agecomp.csv"
## min & max ages
age_min <- 3
age_max <- 6
## years, if any, of age-comp to skip; see below
age_skip <- 0
## 3. file with harvest data
## [n_yrs x 2] matrix of obs catch; 1st col is calendar yr
fn_harv <- "chin_harv.csv"
## 4. covariate metadata
cov_meta_file <- "chin_cov_metadata.csv"
## URL for example data files
## set to NULL if using a local folder/directory
ex_url <- "https://raw.githubusercontent.com/mdscheuerell/willamette/master/data/"Next we specify some additional information regarding scenarios and forecasting, if desired.
## number of years of forecasts
n_fore <- 1
## file where to save JAGS model
fn_jags <- "Willamette_Chin_SR_flow_models_mainstem_JAGS.txt"
## upper threshold for Gelman & Rubin's (1992) potential scale reduction factor (Rhat).
Rhat_thresh <- 1.1
## threshold values for CI's
CI_vec <- c(0.025,0.5,0.975)
## FLOW SCENARIOS
## dir with flow scenarios
scen_dir <- "WBR_ScenarioPreds_06Mar2017"
## scenario names
scen_names <- c("unregulated","BiOp","lowApril","Recess",
"LowRecess","BiOpTrib","LowSummer","NMFS1",
"Corps1","Corps2","OSU1")
## spawner abundances to consider
scen_sp <- c(5,10,15,20)*1000
## first year of scenarios
scen_yr_frst <- 1929
## last year of scenarios
scen_yr_last <- 20073 Function definitions
## Re2prec
## converts any real number to round/floor/ceiling
Re2prec <- function(x,map="round",prec=1) {
## 'map' can be round, floor, or ceiling
## 'prec' is nearest value (eg, 0.1 means to nearest tenth); default 1 gives normal behavior
if(prec<=0) { stop("\"prec\" cannot be less than or equal to 0") }
do.call(map,list(x/prec))*prec
}
## over90k
## returns the total number of days that flow exceeded 90k cfs
over90k <- function(x) {
return(sum(x >= 90000))
}
## min7mean
## returns the min of 7-day means over period
min7mean <- function(x) {
return(round(min(filter(x, rep(1,7)/7, "convolution", sides=1), na.rm=TRUE),0))
}
## max7mean
## returns the max of 7-day means over period
max7mean <- function(x) {
return(round(max(filter(x, rep(1,7)/7, "convolution", sides=1), na.rm=TRUE),0))
}
## med7mean
## returns the max of 7-day means over period
med7mean <- function(x) {
return(round(median(filter(x, rep(1,7)/7, "convolution", sides=1), na.rm=TRUE),0))
}
## rng7mean
## returns the max of 7-day means over period
rng7mean <- function(x) {
return(round(diff(range(filter(x, rep(1,7)/7, "convolution", sides=1), na.rm=TRUE)),0))
}
## range2
## returns the range as a scalar = max-min
range2 <- function(x) {
return(diff(range(x)))
}
## DD15
## returns the cumulative number of degree days > 15C
DD15 <- function(x) {
return(sum(filter(x, rep(1,7)/7, "convolution", sides=1) >= 15, na.rm=TRUE))
}
## colVars; from Gelman
## returns the column-wise variance of a matrix
colVars <- function(a) {
n <- dim(a)[[1]]
c <- dim(a)[[2]]
mm <- matrix(.colMeans(a, n, c), n, c, byrow = TRUE)
return(.colMeans(((a - mm) ^ 2), n, c) * n / (n - 1))
}
## waic; from Gelman
## computes WAIC based on pointwise log-like
waic <- function(log_lik) {
S <- nrow(log_lik)
n <- ncol(log_lik)
lpd <- log(colMeans(exp(log_lik)))
p_waic <- colVars(log_lik)
elpd_waic <- lpd - p_waic
waic <- -2*elpd_waic
loo_weights_raw <- 1/exp(log_lik-max(log_lik))
loo_weights_normalized <- loo_weights_raw /
matrix(colMeans(loo_weights_raw),nrow=S,ncol=n,byrow=TRUE)
loo_weights_regularized <- pmin(loo_weights_normalized, sqrt(S))
elpd_loo <- log(colMeans(exp(log_lik)*loo_weights_regularized) /
colMeans(loo_weights_regularized))
p_loo <- lpd - elpd_loo
pointwise <- cbind(waic,lpd,p_waic,elpd_waic,p_loo,elpd_loo)
total <- colSums(pointwise)
se <- sqrt(n*colVars(pointwise))
return(list(waic=total["waic"],
elpd_waic=total["elpd_waic"],
p_waic=total["p_waic"],
elpd_loo=total["elpd_loo"],
p_loo=total["p_loo"],
pointwise=pointwise,
total=total,
se=se))
}4 Loading the fish data
Here we load in the first three data files and do some simple calculations and manipulations.
First the spawner data:
## escapement
dat_esc <- read.csv(text = getURL(paste0(ex_url,fn_esc)))
## use total counts
dat_esc <- dat_esc[dat_esc$group=="total",-1]
## inspect
dat_esc## year unmarked
## 35 1999 7689
## 36 2000 9703
## 37 2001 11904
## 38 2002 21538
## 39 2003 22767
## 40 2004 19056
## 41 2005 8665
## 42 2006 9810
## 43 2007 7903
## 44 2008 5463
## 45 2009 7706
## 46 2010 14148
## 47 2011 13427
## 48 2012 8715
## 49 2013 7352
## 50 2014 6696
## 51 2015 9429## years of data
dat_yrs <- dat_esc$year
## number of years of data
n_yrs <- length(dat_yrs)
## get first & last years
yr_frst <- min(dat_yrs)
yr_last <- max(dat_yrs)
## log of escapement
ln_dat_esc <- c(log(dat_esc[,-1]),rep(NA,n_fore))Next the age composition data:
## age comp data
dat_age <- read.csv(text = getURL(paste0(ex_url,fn_age)))
## use total age counts by year
dat_age <- dat_age[dat_age$site=="all" & dat_age$year>=yr_frst & dat_age$year<=yr_last,]
## drop first age_min+age_skip rows; drop site & year col
dat_age <- dat_age[-(1:(age_min+age_skip)),-c(1,2)]
## num of age classes
A <- age_max-age_min+1
## add row(s) of NA's for forecast years
dat_age <- rbind(dat_age,matrix(0,n_fore,A,dimnames=list(n_yrs+seq(n_fore),colnames(dat_age))))
## total num of age obs by cal yr
dat_age[,"sum"] <- apply(dat_age,1,sum)
## row indices for any years with no obs age comp
idx_NA_yrs <- which(dat_age$sum<A,TRUE)
## replace 0's in yrs w/o any obs with NA's
dat_age[idx_NA_yrs,(1:A)] <- NA
## change total in yrs w/o any obs from 0 to A to help dmulti()
dat_age[idx_NA_yrs,"sum"] <- A
## convert class
dat_age <- as.matrix(dat_age)
## inspect
dat_age## age_3 age_4 age_5 age_6 sum
## 127 8 397 359 5 769
## 128 2 158 288 13 461
## 129 1 208 377 6 592
## 130 8 82 241 11 342
## 131 7 173 142 4 326
## 132 29 251 567 22 869
## 133 105 488 280 13 886
## 134 12 467 199 18 696
## 135 30 200 219 5 454
## 136 17 814 303 11 1145
## 137 5 90 288 22 405
## 138 57 413 239 12 721
## 139 96 940 447 11 1494
## 140 66 713 176 7 962
## 18 NA NA NA NA 4And then the harvest data:
## harvest
dat_harv <- read.csv(text = getURL(paste0(ex_url,fn_harv)))
## trim to correct years & drop year col
dat_harv <- c(dat_harv[dat_harv$year>=yr_frst & dat_harv$year<=yr_last,-1],rep(0,n_fore))5 Loading the covariates
The analyses are based upon several environmental indicators related to river discharge and temperature. Load the data file containing all of the metadata regarding the specific covariates to be used.
cov_meta <- read.csv(text = getURL(paste0(ex_url,cov_meta_file)), stringsAsFactors = FALSE)
cov_meta$code <- gsub("\"","",cov_meta$code)
cov_meta$begin <- gsub("\"","",cov_meta$begin)
cov_meta$end <- gsub("\"","",cov_meta$end)We need to define the beginning and ending dates for the covariates
yr1 <- yr_frst
yr2 <- yr_last - age_min + n_fore
## start date
startDate <- paste0(yr1 + min(cov_meta$lag_1),"-01-01")
## end date
endDate <- paste0(yr2 + max(cov_meta$lag_2),"-12-31")5.1 River discharge
We begin by getting the daily flow data from the US Geological Service National Water Information System for the complete time period of interest.
## metadata for flow covariates
flow_meta <- subset(cov_meta, covariate=="flow" & flag==1)
## data to get: flow (cfs)
parameterCD <- "00060"
## get all flow data for period of interest & gages
gages <- unique(flow_meta$gage)
n_gages <- length(gages)
tmp <- readNWISDaily(gages[1], parameterCD, startDate, endDate, convert=FALSE)
tmp$yr <- floor(tmp$DecYear)
flow_data <- tmp[,c("Date","waterYear","yr","Month","Day")]
colnames(flow_data) <- c("date","H2Oyr","yr","mon","day")
flow_data[,as.character(gages[1])] <- tmp$Q
if(n_gages > 1) {
for(i in 2:n_gages) {
tmp <- readNWISDaily(gages[i], parameterCD, startDate, endDate, convert=FALSE)
flow_data[,as.character(gages[i])] <- tmp$Q
}
}## There are 6574 data points, and 6574 days.Now we can extract the specific flow covariates that relate to each of the hypotheses about the affected life stage.
cov_flow <- matrix(NA,length(seq(yr1,yr2)),dim(flow_meta)[1]+1)
n_mods <- dim(cov_flow)[2] - 1
cov_flow[,1] <- seq(yr_frst,yr2)
for(i in 1:dim(flow_meta)[1]) {
fn <- get(flow_meta[i,"short_name"])
cnt <- 1
for(t in yr1:yr2) {
beg <- paste0(t+flow_meta[i,"lag_1"],"-",flow_meta[i,"begin"])
end <- paste0(t+flow_meta[i,"lag_2"],"-",flow_meta[i,"end"])
tmp <- flow_data[flow_data$date>=beg & flow_data$date<=end, as.character(flow_meta[i,"gage"])]
cov_flow[cnt,i+1] <- fn(tmp)
cnt <- cnt+1
}
}
print(round(cov_flow,0))## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 1999 117100 53657 15014 21214 36129 21114 19443 22157 36129 7354
## [2,] 2000 81000 36371 9681 19586 30286 20604 18900 17000 30286 6854
## [3,] 2001 26700 11729 6667 14871 20143 13476 14114 14643 18371 5176
## [4,] 2002 81943 35143 8467 15629 54300 45833 17671 15429 16271 6469
## [5,] 2003 74100 29157 7089 16229 39429 32340 34014 14843 21471 6016
## [6,] 2004 75943 31243 8930 15257 19686 10756 15343 14414 18600 6803
## [7,] 2005 54343 12829 9597 17029 27286 17689 16443 15414 27286 7043
## [8,] 2006 125286 26571 9583 17243 25929 16346 18571 15929 19071 6497
## [9,] 2007 83557 40271 7937 15143 19529 11591 18043 14286 15757 5854
## [10,] 2008 73243 34250 16071 24186 43214 27143 17114 24400 43214 6570
## [11,] 2009 95286 22457 8699 19714 39557 30859 18543 17271 39557 6310
## [12,] 2010 61129 21100 13514 24486 71157 57643 20214 18029 29214 6581
## [13,] 2011 84286 37814 16371 28657 49514 33143 31314 22257 29986 7639
## [14,] 2012 122871 29071 15186 26614 93357 78171 33457 19229 41200 7877
## [15,] 2013 89157 24286 10257 15600 29300 19043 16529 14457 19043 7047
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 8490 81000 36371 9681 19586 30286 20604 18900 30286 17000 6854
## [2,] 7334 26700 11729 6667 14871 20143 13476 14114 18371 14643 5176
## [3,] 5659 81943 35143 8467 15629 54300 45833 17671 16271 15429 6469
## [4,] 6673 74100 29157 7089 16229 39429 32340 34014 21471 14843 6016
## [5,] 6246 75943 31243 8930 15257 19686 10756 15343 18600 14414 6803
## [6,] 7225 54343 12829 9597 17029 27286 17689 16443 27286 15414 7043
## [7,] 7327 125286 26571 9583 17243 25929 16346 18571 19071 15929 6497
## [8,] 6809 83557 40271 7937 15143 19529 11591 18043 15757 14286 5854
## [9,] 6375 73243 34250 16071 24186 43214 27143 17114 43214 24400 6570
## [10,] 7602 95286 22457 8699 19714 39557 30859 18543 39557 17271 6310
## [11,] 7203 61129 21100 13514 24486 71157 57643 20214 29214 18029 6581
## [12,] 8609 84286 37814 16371 28657 49514 33143 31314 29986 22257 7639
## [13,] 9646 122871 29071 15186 26614 93357 78171 33457 41200 19229 7877
## [14,] 8733 89157 24286 10257 15600 29300 19043 16529 19043 14457 7047
## [15,] 7601 105286 23829 11100 19514 38543 27443 17557 34843 15329 7426
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## [1,] 7334 10000 3146 20143 12304 8747 11396 20143 17179 14114 6029
## [2,] 5659 6701 1526 54300 25957 15486 38814 54300 29214 17671 36629
## [3,] 6673 7591 1123 67600 36043 14643 52957 39429 35843 34014 5414
## [4,] 6246 7287 1271 63186 19029 13243 49943 19300 16543 15343 3957
## [5,] 7225 11643 4840 35371 9490 6543 28829 22771 21114 16443 6329
## [6,] 7327 9090 2047 87900 21029 14186 73714 25929 22921 18571 7357
## [7,] 6809 8594 2097 44200 23943 11743 32457 19529 18800 18043 1486
## [8,] 6375 7544 1690 49929 22264 16886 33043 24114 21179 17114 7000
## [9,] 7602 13800 7230 32829 23086 10571 22257 26543 22721 18543 8000
## [10,] 7203 9021 2711 55129 19643 13114 42014 48671 25543 20214 28457
## [11,] 8609 11571 4990 49657 35500 13471 36186 49514 42214 31314 18200
## [12,] 9646 13243 5604 98929 40321 19971 78957 93357 44064 33457 59900
## [13,] 8733 12500 4623 31043 18871 15171 15871 29300 19886 16529 12771
## [14,] 7601 17029 9981 105286 41529 17557 87729 38543 25336 17557 20986
## [15,] 7915 9560 2134 55800 14529 10104 45696 14586 13707 10986 3600write.csv(cov_flow, row.names = FALSE,
file="Willamette_Chin_SR_flow_models_mainstem_covariates.csv")Here are plots of the pairwise correlation coefficients for each of the flow covariates.
## min cor to consider
thrsh <- 0.7
## color palette
clr <- viridis(n=11, begin=0.1, end=0.9, option="B")
## transformed matrix of corr coefs
tcc <- t(apply(cor(cov_flow[,-1]),1,rev))
nn <- dim(tcc)[1]
par(mai=rep(0.5,4), omi=c(0,0,0,0.5))
## color image
image(tcc, xaxt="n", yaxt="n", col=clr)
axis(1,at=(seq(0,nn,length.out=nn)+1/nn)/nn, labels=seq(nn), cex.axis=0.7, tick=FALSE, line=-1)
axis(2,at=rev(seq(0,nn,length.out=nn)+1/nn)/nn, labels=seq(nn), las=1, cex.axis=0.7, tick=FALSE, line=-0.8)
axis(3,at=(seq(0,nn,length.out=nn)+1/nn)/nn, labels=seq(nn), cex.axis=0.7, tick=FALSE, line=-1)
axis(4,at=rev(seq(0,nn,length.out=nn)+1/nn)/nn, labels=seq(nn), las=1, cex.axis=0.7, tick=FALSE, line=-0.8)
grid(nx=n_mods,ny=n_mods, lty="solid")
## show major breaks by life stage
bls <- cumsum(table(cov_meta$group))[-length(unique(cov_meta$group))] + 1
abline(h=seq(par()$usr[1],par()$usr[2],length.out=nn+1)[nn-bls+2], lwd=4, col="white")
abline(v=seq(par()$usr[1],par()$usr[2],length.out=nn+1)[bls], lwd=4, col="white")
## add legend
legend(x="right", legend=rev(seq(0,10)/10), title="Corr",
inset=c(-0.15,0), xpd=NA,
col=rev(viridis(n=11, begin=0.1, end=0.9, option="B")), bg="white",
pch=15, pt.cex=1.2, cex=0.7, bty="n")
5.2 Water temperature
Next we get the daily temperature data from the US Geological Service National Water Information System for the complete time period of interest.
## metadata for flow covariates
temp_meta <- subset(cov_meta, covariate=="temp" & flag==1)
## data to get: temp (C)
parameterCD <- "00010"
## get all flow data for period of interest & gages
gages <- unique(temp_meta$gage)
n_gages <- length(gages)
tmp <- readNWISDaily(gages[1], parameterCD, startDate, endDate, convert=FALSE)
tmp$yr <- floor(tmp$DecYear)
temp_data <- tmp[,c("Date","waterYear","yr","Month","Day")]
colnames(temp_data) <- c("date","H2Oyr","yr","mon","day")
temp_data[,as.character(gages[1])] <- tmp$Q
if(n_gages > 1) {
for(i in 2:n_gages) {
tmp <- readNWISDaily(gages[i], parameterCD, startDate, endDate, convert=FALSE)
flow_data[,as.character(gages[i])] <- tmp$Q
}
}Now we can extract the specific temperature covariates that relate to each of the hypotheses about the affected life stage.
cov_temp <- matrix(NA,length(seq(yr1,yr2)),dim(temp_meta)[1]+1)
cov_temp[,1] <- seq(yr_frst,yr2)
for(i in 1:dim(temp_meta)[1]) {
fn <- get(temp_meta[i,"short_name"])
cnt <- 1
for(t in yr1:yr2) {
beg <- paste0(t+temp_meta[i,"lag_1"],"-",temp_meta[i,"begin"])
end <- paste0(t+temp_meta[i,"lag_2"],"-",temp_meta[i,"end"])
tmp <- temp_data[temp_data$date>=beg & temp_data$date<=end, as.character(temp_meta[i,"gage"])]
if(length(tmp)>0) {
cov_temp[cnt,i+1] <- fn(tmp)
} else {
cov_temp[cnt,i+1] <- NA
}
cnt <- cnt+1
}
}
cov_tempUnfortunately, there are missing years of temperature data, even at the Portland gage.
6 Specifying model in JAGS
Now we can specify the model in JAGS. Note that the code below is not written to be universally generic with respect to the number of covariates, but rather to emphasize how to incorporate the three we used in this specific application. The important thing is the number of covariate parameters in the PRIORS and LIKELIHOOD sections (i.e., there must be a unique c_Name parameter for each of the MM covariates).
cat("
model {
##--------
## PRIORS
##--------
## alpha = exp(a) = intrinsic productivity
mu_Rkr_a ~ dnorm(0,1e-3)I(-4,4);
alpha <- exp(mu_Rkr_a);
E_Rkr_a <- mu_Rkr_a + var_Qr/(2 - 2*phi^2);
## covariate effects
c1 ~ dnorm(0,0.001)
c2 ~ dnorm(0,0.001)
## strength of dens depend
Rkr_b ~ dunif(0,0.01);
## AR(1) coef for proc errors
phi ~ dunif(-0.999,0.999);
## process variance for recruits model
sd_Qr ~ dunif(0.001,20);
tau_Qr <- pow(sd_Qr,-2);
var_Qr <- pow(sd_Qr,2)
## innovation in first year
innov_1 ~ dnorm(0,tau_Qr*(1-phi*phi));
## obs variance for spawners
sd_Rs ~ dunif(0.001,3);
tau_Rs <- pow(sd_Rs,-2);
var_Rs <- pow(sd_Rs,2)
## unprojectable early recruits;
## hyper mean across all popns
Rec_mu ~ dnorm(0,0.001);
## hyper SD across all popns
Rec_sig ~ dunif(0,100);
## precision across all popns
Rec_tau <- pow(Rec_sig,-2);
## multipliers for unobservable total runs
ttl_run_mu ~ dunif(1,5);
ttl_run_tau ~ dunif(1,20);
## maturity schedule
## unif vec for Dirch prior
for(i in 1:A) { theta[i] <- 1 }
## hyper-mean for maturity
p_mu ~ ddirch(theta);
## hyper-prec for maturity
p_pi ~ dunif(0.001,1e3);
for(t in 1:(n_yrs-age_min+n_fore)) { p_vec[t,1:A] ~ ddirch(p_mu*p_pi) }
##------------
## LIKELIHOOD
##------------
## 1st brood yr requires different innovation
## predicted recruits in BY t
ln_Rkr_a[1] <- mu_Rkr_a + c1*covar[1] + c2*covar[1]^2;
E_ln_Rec[1] <- ln_Sp[1] + ln_Rkr_a[1] - Rkr_b*Sp[1] + phi*innov_1;
tot_ln_Rec[1] ~ dnorm(E_ln_Rec[1],tau_Qr);
res_ln_Rec[1] <- tot_ln_Rec[1] - E_ln_Rec[1];
## median of total recruits
tot_Rec[1] <- exp(tot_ln_Rec[1]);
## R/S
ln_RS[1] <- tot_ln_Rec[1] - ln_Sp[1];
## brood-yr recruits by age
for(a in 1:A) {
Rec[1,a] <- max(1,tot_Rec[1] * p_vec[1,a]);
}
## brood years 2:(n_yrs-age_min)
for(t in 2:(n_yrs-age_min+n_fore)) {
## predicted recruits in BY t
ln_Rkr_a[t] <- mu_Rkr_a + c1*covar[t] + c2*covar[t]^2;
E_ln_Rec[t] <- ln_Sp[t] + ln_Rkr_a[t] - Rkr_b*Sp[t] + phi*res_ln_Rec[t-1];
tot_ln_Rec[t] ~ dnorm(E_ln_Rec[t],tau_Qr);
res_ln_Rec[t] <- tot_ln_Rec[t] - E_ln_Rec[t];
## median of total recruits
tot_Rec[t] <- exp(tot_ln_Rec[t]);
## R/S
ln_RS[t] <- tot_ln_Rec[t] - ln_Sp[t];
## brood-yr recruits by age
for(a in 1:A) {
Rec[t,a] <- max(1,tot_Rec[t] * p_vec[t,a]);
}
} ## end t loop over year
## get total cal yr returns for first age_min yrs
for(i in 1:(age_min+age_skip)) {
ln_tot_Run[i] ~ dnorm(ttl_run_mu*Rec_mu,Rec_tau/ttl_run_tau);
tot_Run[i] <- exp(ln_tot_Run[i]);
}
## get predicted calendar year returns by age
## matrix Run has dim [(n_yrs-age_min) x A]
## step 1: incomplete early broods
## first cal yr of this grp is first brood yr + age_min + age_skip
for(i in 1:(age_max-age_min-age_skip)) {
## projected recruits
for(a in 1:(i+age_skip)) {
Run[i,a] <- Rec[(age_skip+i)-a+1,a];
}
## imputed recruits
for(a in (i+1+age_skip):A) {
lnRec[i,a] ~ dnorm(Rec_mu,Rec_tau);
Run[i,a] <- exp(lnRec[i,a]);
}
## total run size
tot_Run[i+age_min+age_skip] <- sum(Run[i,1:A]);
## predicted age-prop vec for multinom
for(a in 1:A) {
age_v[i,a] <- Run[i,a] / tot_Run[i+age_min];
}
## multinomial for age comp
dat_age[i,1:A] ~ dmulti(age_v[i,1:A],dat_age[i,A+1]);
lp_age[i] <- logdensity.multi(dat_age[i,1:A],age_v[i,1:A],dat_age[i,A+1]);
}
## step 2: info from complete broods
## first cal yr of this grp is first brood yr + age_max
for(i in (A-age_skip):(n_yrs-age_min-age_skip+n_fore)) {
for(a in 1:A) {
Run[i,a] <- Rec[(age_skip+i)-a+1,a];
}
## total run size
tot_Run[i+age_min+age_skip] <- sum(Run[i,1:A]);
## predicted age-prop vec for multinom
for(a in 1:A) {
age_v[i,a] <- Run[i,a] / tot_Run[i+age_min];
}
## multinomial for age comp
dat_age[i,1:A] ~ dmulti(age_v[i,1:A],dat_age[i,A+1]);
lp_age[i] <- logdensity.multi(dat_age[i,1:A],age_v[i,1:A],dat_age[i,A+1]);
}
## get predicted calendar year spawners
## first cal yr is first brood yr
for(t in 1:(n_yrs+n_fore)) {
## obs model for spawners
Sp[t] <- max(1,tot_Run[t] - dat_harv[t]);
ln_Sp[t] <- log(Sp[t]);
ln_dat_esc[t] ~ dnorm(ln_Sp[t], tau_Rs);
lp_esc[t] <- logdensity.norm(ln_dat_esc[t],ln_Sp[t], tau_Rs)
}
} ## end model description
", file=fn_jags)7 Fitting the model
The last thing we need to do before fitting the model in JAGS is to specify:
- the data and indices that go into the model;
- the model parameters and states that we want JAGS to return;
- the MCMC control parameters.
Please note that the following code takes ~3 min to run on a quad-core machine with 3.5 GHz Intel processors.
## 1. data to pass to JAGS
dat_jags <- c("dat_age","ln_dat_esc","dat_harv","covar",
"n_yrs","A","age_min","age_max","age_skip","n_fore")
## 2. model params/states for JAGS to return
par_jags <- c("alpha","E_Rkr_a","mu_Rkr_a","Rkr_b","ln_Rkr_a",
"c1","c2","lp_age","lp_esc",
"Sp","Rec","tot_ln_Rec","ln_RS","p_vec",
"var_Qr","var_Rs","res_ln_Rec")
## 3. MCMC control params
## MCMC parameters
mcmc_chains <- 4
mcmc_length <- 4.5e5
mcmc_burn <- 2e5
mcmc_thin <- 200
## total number of MCMC samples
mcmc_samp <- (mcmc_length-mcmc_burn)*mcmc_chains/mcmc_thin
## function to create JAGS inits
init_vals <- function() {
list(mu_Rkr_a=1, c1=rnorm(1,0,0.1), c2=rnorm(1,0,0.1),
Rkr_b=1/exp(mean(ln_dat_esc, na.rm=TRUE)),
p_pi=1, p_mu=rep(1,A),
p_vec=matrix(c(0.04,0.50,0.45,0.01),n_yrs-age_min+n_fore,A,byrow=TRUE),
Rec_mu=log(1000),
Rec_sig=0.1,
tot_ln_Rec=rep(log(1000),n_yrs-age_min+n_fore),
sd_Rs=0.1,
innov_1=0,
phi=0.5)
}
## fit the model in JAGS & store results
mod_fits <- vector("list", n_mods)
for(i in 1:n_mods) {
covar <- (cov_flow[,i+1]-mean(cov_flow[,i+1]))/sd(cov_flow[,i+1])
dat_jags <- c("dat_age","ln_dat_esc","dat_harv","covar",
"n_yrs","A","age_min","age_max","age_skip","n_fore")
## list of model info for JAGS
mod_jags <- list(data=dat_jags,
inits=init_vals,
parameters.to.save=par_jags,
model.file=fn_jags,
n.chains=as.integer(mcmc_chains),
n.iter=as.integer(mcmc_length),
n.burnin=as.integer(mcmc_burn),
n.thin=as.integer(mcmc_thin),
DIC=TRUE)
mod_fits[[i]] <- do.call(jags.parallel, mod_jags)
}
save(list=ls(), file="Willamette_Chin_SR_flow_models_mainstem.RData")8 Scenarios for flow effects
Now we can calculate the anticipated effects for each of the flow covariates under all of the future flow scenarios.
## indices for covars with scenarios
dd <- dir(paste0(getwd(),"/",scen_dir))
fs_idx <- seq(nn) %in% unique(na.omit(as.numeric(unlist(strsplit(dd,"[^0-9]+")))))
## data frame for results
mod_scen <- flow_meta$flow_scen[fs_idx]
n_cov <- length(mod_scen)
## years with flow covariates for predictions
yrs_scen <- c(scen_yr_frst:scen_yr_last)
n_yrs_scen <- length(scen_yr_frst:scen_yr_last)
## number of flow scenarios
n_scen <- length(scen_names)
## number of spawner levels to consider
n_sp_scen <- length(scen_sp)
## fn for predicted R/S
RpS <- function(p,spawners,x) {
return(round(exp(p[,1]-p[,2]*spawners+p[,3]*x+p[,4]*x^2),2))
}
## empty table to hold predictions for each combo of covariate, year, spawner
tbl_scens_all <- expand.grid(yr=yrs_scen, scen=scen_names, sp=scen_sp)
tbl_scens_all <- data.frame(tbl_scens_all[,c("yr","sp","scen")], flow=NA)
tbl_scens_mcmc <- matrix(NA,dim(tbl_scens_all)[1],mcmc_samp)
## reset counter
nt <- n_yrs_scen * n_scen * n_sp_scen
#for(i in mod_scen) {
for(i in 1:n_cov) {
## get model S-R params
pp <- mod_fits[[i]]$BUGSoutput$sims.array[,,c("mu_Rkr_a","Rkr_b","c1","c2")]
dim(pp) <- c(mcmc_samp,mcmc_chains)
## Read in flow values for the selected covariate
flows <- read.csv(paste0(scen_dir,"/AllYears_ScenarioPreds_CovariateNum_",i,".csv"))
flw_z <- rep(melt(scale(flows[,-1]))$value,n_sp_scen)
tbl_scens_all[,"flow"] <- rep(melt(flows,id="year")$value,n_sp_scen)
for(j in 1:nt) {
tbl_scens_mcmc[j,] <- RpS(pp,tbl_scens_all[,"sp"],flw_z[j])
}
## print predictions across mcmc
if(i==1) { apd <- FALSE } else { apd <- TRUE }
tbl_scens_tmp <- data.frame(cov=i,tbl_scens_all,tbl_scens_mcmc)
write.table(tbl_scens_tmp, append=apd, sep=",", row.names=FALSE, col.names=!apd,
file=paste0("Willamette_Chin_RS_pred_all.csv"))
## summary table with 2.5%, 50% & 97.5% of posterior R/S
tbl_scens <- data.frame(tbl_scens_tmp[,c(1:5)],lo=NA,med=NA,hi=NA)
tbl_scens[,-c(1:5)] <- round(t(apply(tbl_scens_mcmc,1,quantile,CI_vec)),2)
write.table(tbl_scens, append=apd, sep=",", row.names=FALSE, col.names=!apd,
file=paste0("Willamette_Chin_RS_pred_smry.csv"))
}